class: inverse,left, middle background-image: url(background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 3: Percepción remota, análisis masivo y Google Earth Engine. ### Pre-procesamiento satelital I Gabriel Castro mail gabriel.castro.b@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Octubre 2024</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ - Imagenes comerciales: Planet Scope Skysat. - Combinaciones de bandas. - Calibración radiométrica. (TOA radianza - TOA reflectancia) - Corrección atmosférica. - Filtros de calidad de la información satelital. - Compensación topográfica por iluminación. ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> © Allison Horst ] --- ### Recordemos. <center> <img src="data:image/png;base64,#Signatura_vegetacion.PNG" alt="drawing" width="800" height="480"/> <center> --- ### 1) Imagenes satelitales PlanetScope - Skysat. #### - [**Información general**](https://assets.planet.com/docs/Planet_Combined_Imagery_Product_Specs_letter_screen.pdf) #### - Especificaciones tecnicas. #### - Niveles de procesamiento. #### - Metadatos e interpretación de la información. <div style="text-align: center;"> <img src="data:image/png;base64,#gif_sat.gif" alt="Descripción del GIF" style="width: 500px;"/> </div> --- <center> <img src="data:image/png;base64,#coleccion.PNG" alt="drawing" width="700" height="650"/> <center> --- <center> <img src="data:image/png;base64,#planet_catalogo.PNG" alt="drawing" width="950" height="500"/> <center> [**Especificaciones catálogo Skysat**](https://developers.planet.com/docs/data/skysat/) --- ### 2) Combinaciones de bandas. A partir de las signaturas espectrales es posible entender el comportamiento de los objetos frente a diferentes longitudes de onda. Esto nos permite crear combinaciones de bandas en los canales **RGB** para resaltar elementos especÃficos en superficie. ```r #librerias necesarias ---- library(satellite) library(tidyverse) library(terra) library(raster) #Producto Analytic - ortho. skysat <- rast("Skysat_ViñadelMar/SkySat_08_feb_2024_vina.tif") par(mfrow = c(1, 3)) plotRGB(skysat, r=3, g=2, b=1, stretch = "lin", axes = TRUE) plotRGB(skysat, r=4, g=2, b=1, stretch = "lin", axes = TRUE) plotRGB(skysat, r=3, g=4, b=1, stretch = "lin", axes = TRUE) dev.off() ``` --- ### 2) Combinaciones de bandas (Skysat). <center> <img src="data:image/png;base64,#comb_bandas_skysat.PNG" alt="drawing" width="1000" height="500"/> <center> --- ### 2) Combinaciones de bandas (Landsat). <center> <img src="data:image/png;base64,#comb_bandas.PNG" alt="drawing" width="1000" height="500"/> <center> --- ### 3) Calibración radiométrica. <center> <img src="data:image/png;base64,#calib_radiometrica.PNG" alt="drawing" width="1000" height="500"/> <center> --- ### 3) Calibración radiométrica. En un gran número de sensores (TM, ETM+, OLI, etc.) se establece la conversión a **radianza espectral** en base a una **relación lineal** existente entre el nivel digital bruto y los valores de radianza (unidades fÃsicas). Sin embargo, **siempre es importante** revisar los manuales de usuario de cada producto con el que estemos trabajando. De tal forma que la ecuación que describe la radianza espectral, según Chander et al. (2009) es: <left> <img src="data:image/png;base64,#formula_rad.PNG" alt="drawing" width="1000" height="300"/> <left> --- ### 3) Calibración radiométrica (skysat) <center> <img src="data:image/png;base64,#radio_skysat.PNG" alt="drawing" width="1000" height="500"/> <center> --- ### 3) Calibración radiométrica (skysat) <center> <img src="data:image/png;base64,#Esun_skysat.PNG" alt="drawing" width="800" height="500"/> <center> --- ### 3) Calibración radiométrica (skysat) #### Niveles digitales a radianza (TOA). ```r # 5.- Calibración radiometrica ND TOA radianza (manual) ---- TOA_rad_skysat <- skysat * 0.01 # factor escalado # Comparar histogramas: par(mfrow = c(1, 2)) hist(skysat[[1]], main = "histograma banda azul ND") hist(TOA_rad_skysat[[1]], main = "histograma banda azul TOA") dev.off() ``` <center> <img src="data:image/png;base64,#histogramas.PNG" alt="drawing" width="900" height="300"/> <center> --- <center> <img src="data:image/png;base64,#ND_TOA_REF.PNG" alt="drawing" width="1050" height="600"/> <center> --- ### 3) Calibración radiométrica #### Niveles digitales a reflectancia aparente (TOA reflectancia) - Transformar la radianza a reflectancia al tope de la atmosfera es un paso intermedio a la corrección atmosférica final. - Si bien la mayorÃa de los métodos de corrección transforman automáticamente los niveles digitales a reflectancia superficial, es importante conocer el **porcentaje** de energÃa que es reflejada por los objetos en la superficie y que finalmente es captada por el sensor. - Al ser la reflectancia al tope de la atmosfera, se le conoce como **reflectancia aparente**. --- ### 3) Calibración radiométrica #### Niveles digitales a reflectancia aparente (TOA reflectancia) <center> <img src="data:image/png;base64,#refp.PNG" alt="drawing" width="1000" height="480"/> <center> --- ### 3) Calibración radiométrica #### Niveles digitales a reflectancia aparente (TOA reflectancia). ```r # 6.- Calibración radiometrica ND TOA reflectancia ---- ## proceso manual -------------------------------------- d <- calcEarthSunDist(date = "2024-02-08", formula = "ESA") # AU ESUN <- c(2009.28, 1820.25, 1583.3, 1114.22) # (W * m2 * sr1 * um1) pi <- pi # 3.14 ## transformar angulos a radianes ---------------------- sun_elev <- 65.4*pi/180 # Sun elev radianes sun_azimuth <- 312.9*pi/180 # Sun azimuth radianes sun_zenith <- (90 - 65.4)*pi/180 # solar zenith TOA_reflectance_skysat <- (pi*TOA_rad_skysat*d^2)/(ESUN*cos(90-sun_elev)) ``` --- <center> <img src="data:image/png;base64,#comparacion.PNG" alt="drawing" width="800" height="600"/> <center> --- ### 4) Corrección atmosferica. <b> Clasificación de correcciones atmosféricas: </b> - **Mediciones <i>in situ</i>** - **Mixtas** Ejemplo: Landsat Ecosystem Disturbance Adaptative Processing System (LEPADS) y Land Surface Reflectance Code (LaSRC) - **A partir de modelos fÃsicos de transferencia radiativa** Ejemplos: Second Simulation of a Satellite Signal in the Solar Spectrum (6S vector version), Multi-angle implementation of atmospheri correction (MAIAC). - **Con datos de la propia imagen** Ejemplos: Simple Dark object Substraction (sDOS), Sen2cor, ACOLITE. --- <center> <img src="data:image/png;base64,#reF_app_SR.PNG" alt="drawing" width="1000" height="600"/> <center> --- ### 4) Corrección atmosferica. <center> <img src="data:image/png;base64,#irradiacion.PNG" alt="drawing" width="700" height="500"/> <center> --- ### 4) Corrección atmosferica. <center> <img src="data:image/png;base64,#softwares.PNG" alt="drawing" width="900" height="500"/> <center> --- ### 4) Corrección atmosferica. <b> Métodos darkest object substraction: </b> En la literatura podemos encontrar diferentes métodos DOS los cuales se basan en principios similares. El paradigma inicial señala que algunos pixeles son cuerpos oscuros cuya reflectancia seria cercana a 0. Esto indica que la radianza captada por el sensor de estos pixeles oscuros, estarÃa influida por el efecto de la **bruma**. Para estos métodos se estima que los piexeles oscuros reflejan al menos entre el 1-2% de la energÃa. (0.01) El método **darkest pixel subtraction** busca reducir o eliminar este efecto de la bruma para obtener la reflectancia superficial. --- ### 4) Corrección atmosferica. <center> <img src="data:image/png;base64,#solarAngles.PNG" alt="drawing" width="920" height="510"/> <center> --- ### 4) Corrección atmosferica. <center> <img src="data:image/png;base64,#DIFF_METODOS2.PNG" alt="drawing" width="960" height="560"/> <center> --- <center> <img src="data:image/png;base64,#METODOS.PNG" alt="drawing" width="730" height="650"/> <center> --- ### Metodo DOS2 Chavez 1996. Este método de corrección por DOS presento una mejora respecto al **<i>simple darkest object subtraction</i>** al tener en consideración la corrección de la transmitancia atmosférica desde la superficie terrestre hacia el sol. También conocido como método CosTz, esta mejora utiliza el ángulo zenital del sol el cual está asociado a la transmitancia atmosférica para las fechas especificas en las cuales se tomó la imagen. Sin embargo este método sigue considerando la Transmitancia de la energÃa del sol a la tierra (Tv) como constante (no pierde energÃa en el viaje de entrada se mantiene con valor 1 como el caso del sDOS). --- ### Comparemos ambos metodos: sDOS 1989 v/s DOS2 1996. #### Parametros iniciales sDOS: ```r ## Corrección atmosferica metodo sDOS ---------- Tz <- 1 Tv <- 1 Edown <- 0 ## ---------------------------------------------- skysat_raster <- skysat %>% stack() B <- calcDODN(skysat_raster[[1]]) G <- calcDODN(skysat_raster[[2]]) R <- calcDODN(skysat_raster[[3]]) NIR <- calcDODN(skysat_raster[[4]]) DNm <- c(B,G,R,NIR) hazeValues <- DNm*0.01 ``` **Funcion calcDODN del paquete satellite calcula de forma automatica los valores de radiancia minimo para cada banda considerando n cantidad de piexeles (20% aproximadamente)** --- ### Calculo parametros finales sDOS 1989: ##### 1) L_001: Calculo de la reflectancia minima del 1%. ##### 2) SHV_rad: valores de bruma para cada banda considerando la influencia atmosferica. ##### 3) L_haze: calculo de la influencia atmosferica en la radianza del sensor. ##### 4) SR_sDOS: Calculo sDOS. ```r L_001 <- 0.01 * (ESUN[1:4] * cos(sun_zenith) * Tz + Edown) * Tv / (pi * d^2) SHV_rad <- hazeValues L_haze <- SHV_rad[1:4] - L_001[1:4] SR_sDOS <- pi * (TOA_rad_skysat[[1:4]] - L_haze[1:4]) * d^2 / (Tv * ESUN[1:4] * cos(sun_zenith) * TAUz) + Edown plotRGB(SR_sDOS, r=3, g=2, b=1, stretch = "lin", axes = TRUE, main = "Reflectancia superficial metodo sDOS") writeRaster(SR_sDOS, filename = "Reflectancia_superficial_metodo_sDOS.tif", overwrite = T) ``` --- ### Calculo parametros finales DOS2 1996: #### 1) Los parametros son los mismos donde lo unico que cambia es Tz. ```r ## Corrección atmosferica metodo sDOS ---------- Tz <- cos(sun_zenith) Tv <- 1 Edown <- 0 ## ---------------------------------------------- L_001 <- 0.01 * (ESUN[1:4] * cos(sun_zenith) * Tz + Edown) * Tv / (pi * d^2) SHV_rad <- hazeValues L_haze <- SHV_rad[1:4] - L_001[1:4] SR_DOS2 <- pi * (TOA_rad_skysat[[1:4]] - L_haze[1:4]) * d^2 / (Tv * ESUN[1:4] * cos(sun_zenith) * TAUz) + Edown plotRGB(SR_DOS2, r=3, g=2, b=1, stretch = "lin", axes = TRUE, main = "Rflectancia superficial metodo DOS2") writeRaster(SR_DOS2, filename = "Reflectancia_superficial_metodo_DOS2.tif", overwrite = T) ``` --- ### Comparativa. ```r ### AED resultados obtenidos par(mfrow = c(2,3)) boxplot(TOA_reflectance_skysat[[2]], outline = F, main = "Boxplot TOA_ref", col = "green", ylim = c(0, 0.30)) boxplot(SR_sDOS[[2]], outline = F, main = "Boxplot sDOS", col = "lightblue", ylim = c(0, 0.30)) boxplot(SR_DOS2[[2]], outline = F, main = "Boxplot DOS2", col = "red", ylim = c(0, 0.30)) # ---------------------------------------------------- hist(TOA_reflectance_skysat[[2]]) hist(SR_sDOS[[2]]) hist(SR_DOS2[[2]]) dev.off() ``` --- ### Comparativa. <center> <img src="data:image/png;base64,#comparacion_metodos.PNG" alt="drawing" width="1000" height="550"/> <center> --- ### Resumén calibración radiometrica y corrección atmosférica. - Ya sea estemos trabajando con colecciones satelitales Landsat, Sentinel, Planet, Worldview etc. **<i>Siempre se deben revisar los manuales de uso de cada producto</i>**. Si al momento de recibir las imagenes estas no cuentan con un manual de uso, es obligatorio exigirlos. - **<i>Siempre revisar los metadatos</i>** de los productos. Es importante tener en cuenta cuales son los parámetros que tenemos y cuales deben ser estimados por otros medios. Esto permitirá definir la lÃnea de procesamiento a seguir. - Siempre evaluar que los resultados obtenidos tengan **<i>coherencia</i>**. Valores sobre o bajo el rango esperado, problemas en la visualización, pueden obedecer a errores en la lÃnea de procesamiento. **<i>Es importante realizar un pequeño AED</i>** para garantizar la integridad de los datos y eliminar artefactos. --- ### BibliografÃa. Thomas Nauss, F.D. (2021) satellite: Handling and Manipulating Remote Sensing Data, Satellite - classes and methods for satellite data. Chavez, Jr, Pat. (1996). Image-Based Atmospheric Corrections - Revisited and Improved. Photogrammetric Engineering and Remote Sensing. 62. 1025-1036. Chavez & P.S, 2015. Image-Based Atmospheric Corrections - Revisited and Improved. , no. January 1996. Moran, M.S., Jackson, R.D., Slater, P.N. y Teiuet, P.M., 1992. Evaluation of Simplified Procedures for Retrieval of Land Surface Reflectance Factors from Satellite Sensor Output. , vol. 184, pp. 169-184. Song, C., Woodcock, C.E., Seto, K.C., Lenney, M.P. y Macomber, S.A., 2000. Classification and Change Detection Using Landsat TM Data : When and How to Correct Atmospheric Effects ? , vol. 4257, no. 00. Wang, Y. , Wang, X. y He, H., 2019. An Improved Dark Object Subtraction Method for Atmospheric Correction of Remote Sensing Images [en lÃnea]. S.l.: Springer Singapore. ISBN 9789811399176. Disponible en: http://dx.doi.org/10.1007/978-981-13-9917-6_41. --- class: inverse middle 